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HEAT  CONDUCTION  IN  THREE  DIMENSIONS 


INTRODUCTION 

There  is  a need  for  accurate  three-dimensional  heat  transport  prediction  in  many  areas 
of  applied  research,  including  solidification  of  molten  materials,  ablation  of  missile  skins 
under  aerodynamic  heating,  vulnerability  of  targets  under  thermal  loads,  and  characteriza- 
tion of  high-power  mirror  surfaces  under  extreme  thermal  loads.  In  some  cases  the  problem 
can  adequately  be  described  by  unidimensional  heat  flux,  but  many  problems  cannot  be 
solved  in  a convenient  and  general  form  without  the  added  dimensions  of  transverse  heat 
flow.  The  ongoing  laser  damage  and  vulnerability  program  in  the  Department  of  Defense 
could  more  effectively  make  use  of  limited  test  schedules  with  the  aid  of  a three-dimensional 
code  combining  accuracy  and  computational  ease.  Despite  the  importance  of  the  topic,  the 
literature  reveals  few  solutions  that  may  be  extended  to  practical  problems  of  the  type  pro- 
posed here. 

In  an  effort  to  satisfy  this  need,  and  in  particular  to  characterize  some  high-power 
water-cooled  mirrors  for  future  use  at  the  Optical  Radiation  Laboratory  at  Chesapeake  Beach, 
Md.,  a standard  finite-difference  approximation  has  been  initiated  to  solve  the  three-dimen- 
sional heat  diffusion  equation 

bT 

— (x,  y,  z,  t)  = «v2T(jc,  y,  z,  t).  (1) 

The  complete  difference  form  of  this  equation  is  used  to  provide  maximum  execution  speed 
for  a modest  lattice  point  density.  Although  there  are  more  elaborate  digital  techniques 
which  would  afford  improved  accuracy  with  the  same  lattice  density,  the  additional  pro- 
gram complexity  does  not  justify  the  slight  increase  in  accuracy  near  the  fusion  front  (see 
for  example  Murray  and  Landis  [1]). 

Comparisons  have  been  made  in  the  limiting  case  of  a “flood-loaded”  sample  with  uni- 
directional heat  flow  to  an  exact  premelting  solution  for  temperature  rise  in  a semi-infinite 
slab  in  one  dimension  and  to  internal  temperature  profiles  derived  from  a highly  refined 
fifty-subdivision  one-dimensional  finite-difference  code  by  S.T.  Hanley  [2].  The  agreement 
is  better  than  10%  both  in  internal  temperature  profile  and  in  bum-through  prediction. 


FINITE-DIFFERENCE  MODEL 

The  material  to  be  heated  or  cooled  is  modeled  to  a rectangular  slab  with  arbitrary 
length,  width,  and  thickness  and  with  the  thermal  load  applied  to  the  front  surface.  Before 
heat  can  diffuse  into  the  interior  of  the  slab,  a thermal  “step”  must  take  place  at  the  front 
surface;  this  is  described  by  the  basic  heat  conduction  equation, 


-XvT  = PqA 

Manuscript  submitted  September  9,  1976. 
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in  which  f0  is  the  incident  heat  flux  and  A is  the  absorption  coefficient  or  fraction  of  inci- 
dent power  density  coupled  into  the  material.  K is  the  thermal  conductivity  for  the  material 
and  is  treated  as  a temperature-dependent  parameter.  The  finite-difference  approximation 
used  for  this  equation  is  given  in  the  forward  difference  form  by 


T(i,  j,  fc,  8 +1) 


p0 

■■  T(i,  j,  ft,  8)  + — AcosOAZ 

A 


where  a cosine  function  has  been  introduced  to  allow  arbitrary  incidence  angles  0 from 
normal  to  the  front  surface.  The  X-Y  plane  is  chosen  to  be  the  front  surface,  with  the 
origin  at  one  comer.  The  subscripts  i,  j,  and  ft  refer  to  lattice  points  along  the  X,  Y,  and 
Z axes,  and  8 refers  to  integral  multiples  of  the  step  At  along  the  time  axis.  At  present 
no  loss  terms  are  explicitly  included  for  heat  carried  away  by  convective  cooling  or  radia- 
tion from  the  front  surface.  Under  many  applications  these  terms  represent  a negligible 
fraction  of  the  heat  input  (<  30  W /cm2  for  front  surface  temperatures  less  than  2000K). 
Loss  terms  of  this  type  are  implicitly  treated  in  temperature-dependent  absorption  coef- 
ficient A,  which  also  accounts  for  changes  in  surface  reflectivity. 

The  complete  difference  approximation  for  the  diffusion  of  heat  into  the  interior  of 
the  slab  in  given  by 

7U  7,  ft,  8+1)  = T(i,  j,  ft,  8)  + — [T(/+l,  /,  ft,  8)  + T(i-1,  j,  ft,  8)] 

AX2 


+ — [3X4  7+1,  ft,  8)  + T(i,  7-1,  ft,  8)] 
VY2 


+ — [3X4  7,  ft+1,  8)  + T(i,  7,  fc-1,  8)] 
A Zz 


- 2a At  + + 

\AX2  AY2  AZ2 


T(i,  7,  ft,  8). 
Z2J 


The  usual  simplification  of  choosing  At  such  that  the  T(i,  j,  ft,  8)  term  vanishes  cannot  be 
made  because  of  the  temperature  dependence  of  the  diffusivity 

<6> 

where  p and  C are  the  density  and  specific  heat,  respectively.  The  code  does  select  a time 
step  At  as  a function  of  AX,  AY,  AZ,  and  the  materials  thermal  parameters  which  will 
provide  a stable  convergent  solution.  Stability  is  taken  in  the  sense  that  short-wave  dis- 
turbances are  damped  out.  A slight  improvement  in  accuracy  in  the  internal  approximate 
solution  to  the  diffusion  equation  was  obtained  by  going  to  smaller  time  steps  than  that 
required  for  stability,  with  the  lattice  density  remaining  fixed  in  a manner  analogous  to 
the  one-dimensional  example  given  by  Ozisik  [31 


j. 
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The  edge  and  back  surfaces  are  independently  constrained  to  satisfy  one  of  two  bound- 
ary conditions.  The  first  is  that  of  no  heat  flow  across  the  boundary,  resulting  in 

- K V T = 0 (6) 

for  all  lattice  points  on  that  surface.  A simple  forward  difference  approximation  to  the 
gradient  is  used  here,  as  in  Eq.  (3).  The  alternative  boundary  condition  is  to  heat  sink  all 
lattice  points  in  that  surface  to  the  ambient  or  initial  slab  temperature. 


At  the  instant  any  lattice  point  reaches  the  melting  temperature  for  the  material,  an 
additional  boundary  condition  is  applied: 


dm 

at 


n n 


(7) 


where  L is  the  latent  heat  of  fusion  for  the  material.  This  equation  is  used  to  follow  the 
position  of  the  liquid-solid  interface  along  the  permissible  axis  directions  ± i,  ± and  ± k. 
Once  a lattice  point  reaches  the  melt  temperature  it  is  pinned  at  the  melt  temperature  until 
the  heat  flow  into  the  unit  cell  about  the  lattice  point  minus  the  heat  flow  out,  as  evaluated 
with  Eq.  (7)  along  the  six  unit  vector  directions,  is  equal  to  the  latent  heat  of  melting  for 
the  mass  of  the  unit  cell.  After  this  balance  has  been  met,  the  lattice  point  temperature  is 
allowed  to  rise.  In  this  manner  the  melt  surface  progresses  through  the  interior  of  the  slab 
as  dictated  by  the  local  thermal  gradients.  This  model  is  quite  powerful  in  dealing  with 
irregular  thermal  load  profiles  on  the  front  surface.  Care  must  be  taken  to  ensure  that  the 
time  steps  are  small  compared  with  the  time  required  to  melt  a unit  cell,  in  order  to  avoid 
excessive  errors  in  fusion  front  travel.  For  the  thermal  loads  considered  here,  this  consider- 
ation is  taken  care  of  by  the  internal  time  step  selection  for  stability. 


A brief  description  of  some  of  the  program  features  will  be  given  here,  and  a more 
detailed  discussion  of  program  construction  and  performance  is  given  in  Appendix  A.  The 
program  has  an  internal  library  of  materials  containing  thermal  parameters  with  temperature 
dependence  and  is  expandable  to  any  material  for  which  the  constants  are  known.  The  lat- 
tice is  made  up  of  10  equally  spaced  points  per  axis  along  the  x and  y axes  and  11  equally 
spaced  points  along  the  z axis,  providing  an  accuracy  of  better  than  10%  (as  will  be  shown 
in  the  next  section).  The  front  surface  thermal  load  is  broken  into  three  selectable  cate- 
gories. A “top  hat”  or  uniform  intensity  distribution  of  arbitrary  diameter  and  total  power 
allows  simulation  of  saturated  laser  amplifier  configurations.  A Gaussian  distribution  trun- 
cated at  an  arbitrary  radius  for  the  1/e2  intensity  allows  modeling  of  laser  TEM^  beam  pro- 
files. The  third  distribution  allows  the  user  to  specify  point  by  point  an  arbitrary  thermal 
load  application.  This  flexibility  of  input  distribution  allows  modeling  of  a wide  variety  of 
problems  with  a minimum  of  setup  time. 


APPLICATIONS 

One  of  the  first  uses  of  the  code  was  to  scale  high-reflectance  water-cooled  molybdenum 
mirrors  with  thermal  load.  An  example  of  the  predicted  results  is  shown  in  Fig.  1.  In  the 
model,  a truncated  Gaussian  beam  of  4-cm  radius  impinges  on  a molybdenum  slab  of  0.05-cm 
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TIME  In  MLLISJCOHDS 

Fig.  1 — Thermal  response  of  water-cooled  molybdenum  mirror.  The  front  surface  thermal  re- 
sponse vs  time  is  shown  for  a 98%  reflective  water-cooled  molybdenum  mirror  at  two  power 
levels.  The  thermal  load  is  produced  by  an  8-cm-diam  Gaussian  beam  truncated  at  the  1/e2 
intensity  radius  and  having  100-kW  and  1000-kW  total  power  to  produce  the  upper  and  lower 
curves,  respectively.  Faceplate  thickness  is  0.05  cm. 
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thickness  with  the  back  surface  and  edges  heat  sinked  to  ambient  temperature.  The  front 
surface  reflectivity  is  set  at  a fixed  98%  to  allow  for  aging  and  accumulation  of  impurities. 
Although  enhanced  silver  coatings  are  better  than  this  [4],  it  is  difficult  to  maintain  higher 
reflectance  than  this  over  a period  of  time.  Two  power  levels  are  shown;  the  lower  one 
represents  nominal  power  densities  commonly  in  use  on  high-power  water-cooled  mirrors, 
and  the  upper  curve  represents  an  extrapolation  to  potentially  attainable  power  densities. 
The  cooling  passages  in  the  mirror  wall  allow  some  water  temperature  rise  due  to  the  finite 
flow  rate  through  the  mirror.  At  the  higher  power  level,  the  water  temperature  rise  from 
inlet  to  outlet  may  reach  10°C.  This  effect  has  not  been  included  in  the  model,  however, 
since  it  depends  on  the  design  of  flow  passages.  The  0.05-cm  face  plate  thickness  repre- 
sents that  found  in  some  of  the  better  chemically  etched  mirrors  on  the  market  today.  The 
time  required  to  reach  stable  thermal  conditions  is  short  at  both  power  levels.  Figure  2 
shows  an  overhead  view  of  the  molybdenum  slab  of  Fig.  1,  with  the  temperature  axis  nor- 
mal to  the  front  surface.  The  incident  beam  power  differs  slightly  from  1000  kW  due  to 
the  finite  number  of  lattice  sites  within  the  beam  diameter.  Internal  program  algorithms 
were  designed  to  preserve  flux  densities  specified  by  total  power  and  beam  diameter  and 
then  print  out  the  quantized  beam  power,  rather  than  scale  flux  density  to  total  power. 

A check  of  premelting  accuracy  was  made  by  flood  loading  a slab  of  aluminum  with 
a uniform  beam  distribution  for  times  shorter  than  that  required  for  heat  diffusion  to  the 
back  surface  and  comparing  the  internal  temperature  profile  to  the  exact  one-dimensional 
analytical  result 


T(z,  t)  = 70  + 2 AP0 


(8) 


for  a semi-infinte  slab  of  polished  aluminum.  The  results  after  0.05,  0.1,  0.15,  and  0.2  s 
of  irradiation  are  shown  in  Fig.  3.  For  purposes  of  comparison  the  temperature  depen- 
dence of  the  thermal  parameters  in  the  three-dimensional  code  was  suppressed.  The  upper 
curve  of  each  pair  is  the  three-dimensional  finite-difference  approximation  to  the  exact 
one-dimensional  analytical  solution  given  in  the  lower  curve.  For  all  depths  into  the  slab 
the  accuracy  is  better  than  10%.  Several  alternatives  for  improving  accuracy  are  discussed 
in  the  subsequent  section. 


A limitation  of  thermal  profiling  models  in  the  literature  that  claims  high  accuracy, 
such  as  the  one-dimensional  model  of  Fuhs  and  Fuhs  [5],  is  the  assumption  of  constant 
thermal  parameters.  To  show  the  nature  of  this  restriction,  we  irradiated  a slab  of  polished 
aluminum  with  20,000  W/cm2  for  0.3  and  1.5  s.  First  with  the  temperature  dependence 
contained  in  the  thermal  parameters  and  then  with  thermal  dependence  suppressed.  The 
upper  curve  of  each  pair  in  Fig.  4 shows  the  error  encountered  by  neglecting  the  tempera- 
ture dependence.  As  might  be  expected,  the  error  grows  with  increasing  temperature  and 
becomes  appreciable  after  a modest  amount  of  heating.  In  work  on  laser  damage  of  mate- 
rials such  as  the  AVCO  Everett  study  [6]  done  for  the  Air  Force  Materials  Laboratory, 
significant  improvement  in  modeling  of  experimental  data  should  be  obtained  by  including 
the  effects  of  nonconstant  specific  heat,  thermal  conductivity,  density,  and  absorption. 
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FRONT  SURFACE  CENTIMETERS 

Fig.  2 — Surface  plot  of  water-cooled  molybdenum  mirror  under  thermal  load.  Front  surface  temperatures, 
for  a water-cooled  molybdenum  mirror  under  the  conditions  of  Fig.  1 and  at  the  1000-kW  power  level. 
Shown  after  10  ms  of  radition. 
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DEPTH  FROM  FPONT  SURFACE  IN  CENTI NETEFS 

Fig.  3 — Comparison  of  three-dimensional  code  and  one-dimensional  analytical  solution.  The 
upper  curve  of  each  pair  gives  internal  temperatures  for  a 1-cm-thick  slab  of  aluminum  2024 
uniformly  flood  loaded  at  20000  W/cm2  incident  flux  after  0.05,  0.1,  0.15,  and  0.2  s of  irradia- 
tion. Temperature  dependence  in  all  thermal  parameters  of  the  three-dimensional  code  are  sup- 
pressed. The  lower  curve  of  each  pair  represents  the  exact  one-dimensional  analytical  premelting 
solution  for  the  same  power  density  and  at  the  respective  times  of  irradiation. 
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PERCENT  DEPTH  INTO  SLAB  PROM  FRONT  SURFACE 


Fig.  4 — Temperature  dependence  in  thermal  parameters.  Internal  temperatures  are  shown  for  a 
1-cm-thick  slab  of  polished  aluminum  2024  irradiated  with  20000  W/cm2  in  a flood-loaded  uni- 
form distribution  for  0.3  and  1.5  s.  The  temperature  dependence  of  thermal  parameters  was 
suppressed  in  the  upper  curve  of  each  pair  to  show  the  error  introduced  by  neglecting  tempera- 
ture dependence. 
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Experimental  measurements  of  burn-through  times  on  aluminum  panels  painted  with 
special  “nonleafing”  paint  have  been  made  by  Dr.  R.  Wenzel  of  the  Naval  Research  Lab- 
oratory (NRL).  The  samples  were  mounted  vertically  and  irradiated  at  near  normal  incidence. 
Video  coverage  revealed  retention  of  some  liquid  on  the  surface  throughout  the  melting 
process.  The  measurements  yielded  a statistically  averaged  bum-through  constant  for 
given  target  thickness  and  irradiating  peak  power  density.  R.B.  Brown,  of  NRL,  performed 
absorption  measurements  on  charred  nonleafing  paint  samples  and  rough-melted  aluminum 
samples  at  10.6-//m  wavelengths.  To  test  all  of  the  features  of  the  three-dimensional  code, 
we  ran  a simulation  of  the  painted  aluminum  sample.  Using  for  inputs  to  the  three- 
dimensional  temperature  profiler  Brown’s  measured  absorption  values  for  premelting  and 
postmelting,  temperature-dependent  thermal  parameters  for  aluminum  2024,  a slab  thick- 
ness and  power  density  equal  to  that  used  by  Wenzel  in  his  measurements,  and  45%  of 
total  slab  thickness  retained  as  liquid,  we  heated  the  sample  from  20°C  up  to  melting, 
und  the  melt  interface  progressed  to  the  back  surface  within  5%  of  the  average  burn-through 
time  measured  by  Wenzel. 

As  an  additional  check  on  the  validity  of  the  three-dimensional  code  for  heating  be- 
yond the  onset  of  melting,  a comparison  was  made  to  the  exact  one-dimensional  bum- 
through  expression  in  the  limiting  ease  of  ablation  (complete  melt  removal)  and  constant 
thermal  parameters  derived  by  J.  Rogerson  [7].  Burn  through  occurred  in  5.78  s for  a 
flood-loaded  (beam  diameter  larger  than  sample)  uniform  beam  of  20,000  kW/cm2  inci- 
dent on  a polished  aluminum  2024  slab  1 cm  thick  with  10%  liquid  layer  retention  (mini- 
mum layer  of  3-D  code  at  present)  and  with  temperature  dependence  suppressed.  This 
compared  favorably  with  Rogersons  exact  value  of  5.90  s under  the  same  conditions. 


DISCUSSION 

One  of  the  problems  facing  experimentalists  working  in  the  area  of  laser  interactions 
with  materials  is  the  higher  cost  per  shot  and  longer  down  time  that  are  inherent  in  the 
trend  toward  large  chemical  lasers  as  compared  with  gas  dynamic  lasers.  With  the  increased 
premium  placed  on  experimental  runs,  the  ability  of  this  code  to  predict  trends  and  aid 
selection  of  experiment  parameters  for  maximum  effectiveness  of  limited  data  becomes  a 
valuable  asset. 

An  area  particularly  well  suited  to  modeling  of  the  type  discussed  in  this  paper  is  in 
the  development  of  hot  spot  tracking.  With  the  addition  of  radiative  and  convective  loss 
terms  at  the  surface,  modeling  to  thermal  signatures  from  targets  irradiated  with  arbitrary 
beam  profiles  could  be  performed  using  the  existing  three-dimensional  code.  By  further 
modifying  the  code  to  exclude  certain  of  the  lattice  points  in  the  slab,  geometries  other 
than  flat  surfaces  may  be  used  for  target  simulation.  The  simulated  thermal  signature 
could  be  used  to  evaluate  proposed  hot  spot  tracking  techniques  and  as  an  input  toward 
improving  working  systems. 

At  present  there  are  no  solutions  in  the  literature  capable  of  predicting  spatial  distor- 
tions of  optical  surfaces  under  thermal  loads,  with  the  exception  of  simplified  one- 
dimensional analyses  such  as  that  given  by  Puhs  and  Fuhs  [5].  The  three-dimensional 
code  presented  here  can  in  its  present  state  give  temperature-vs-time  histories  for  arbitrary 
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thermal  loading  and  boundary  conditions,  as  shown  in  Figs.  1 and  2.  With  additions  to 
the  code,  spatial  distortion  in  the  slab  could  be  calculated  directly  from  the  slab  tempera- 
ture matrix  at  any  given  time.  Such  information  could  be  used  in  all  high-power  optics 
configurations. 


The  error  shown  in  Fig.  3,  due  to  10  subdivisions  per  axis  and  a standard  finite- 
difference  approach,  can  be  reduced  in  several  ways.  The  most  straightforward  would 
be  to  increase  the  lattice  point  density  and  thereby  derive  a better  approximation  to  the 
derivatives  in  Eq.  (1).  This  alternative  also  gives  more  surface  lattice  points  with  which 
to  approximate  beam  profiles.  On  the  other  hand  it  would  require  more  computation 
time  and  more  memory  space  from  the  computer.  A second  alternative,  which  was  demon- 
strated in  one-dimension  for  finite-difference  techniques  by  S.T.  Hanley  [2],  is  to  expand 
the  solution  in  spatial  and  time  coordinates  for  improved  accuracy  with  the  same  lattice 
density.  There  will  be  an  increase  in  computation  time  with  this  alternative  due  to  the 
increase  in  mathematical  complexity.  In  the  one-dimensional  case  mentioned  above  it 
was  shown  that  comparable  accuracy  was  achieved  by  either  doubling  the  lattice  density 
or  using  the  coordinate  expansion.  As  shown  in  Fig.  3,  the  nature  of  the  error  is  a small 
offset  or  starting  error  in  the  time  variable;  the  error  remains  essentially  constant  as  heat- 
ing progresses.  An  alternative  procedure  to  increase  accuracy  would  be  to  determine  the 
offset  as  a function  of  lattice  spacing,  time  step,  and  material  parameters  and  start  the 
time  variable  at  a suitable  positive  offset  value.  To  demonstrate  the  improvement  in 
accuracy  obtained  by  this  last  method,  we  determined  the  fixed  time  offset  for  the  exam- 
ple used  in  Fig.  3;  by  starting  the  time  variable  at  this  positive  value  are  obtained  the 
agreement  shown  in  Fig.  5.  The  accuracy  is  better  than  2%  throughout  the  depth  of  the 
slab  and  for  all  times  greater  than  several  offset  steps. 

A particular  experimental  problem  in  the  laser  heating  of  materials  is  measurement 
of  the  liquid  layer  retained  on  the  front  surface  of  the  irradiated  sample  as  bum-through 
progresses.  This  is  difficult  to  measure  directly  and  has  a significant  effect  on  burn- 
through  times.  Figure  6 shows  back  surface  temperature  rise  vs  time  for  a slab  of  alumi- 
num 1 cm  thick  irradiated  with  a 44-kW  Gaussian  beam  of  3-cm  diameter  and  with  10%, 
20%,  40%,  60%,  and  100%  liquid  layers  retained  at  the  front  surface  as  the  melt  pro- 
gresses. The  faster  temperature  rise  corresponds  with  the  thinner  liquid  layer  retained. 

A measurement  could  be  made  with  conventional  back  surface  thermocouple  techniques, 
and  the  melt  layer  could  then  be  derived  by  comparison  to  a chart  such  as  Fig.  6 for 
the  slab  under  study. 
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Fig.  6 — Back  surface  temperature  signature  vs  melt  removal.  The  center  back  surface  tempera- 
ture rise  is  given  for  a 1-cm-thick  painted  aluminum  panel  irradiated  with  a 3 cm  truncated  Gaus- 
sian beam  of  44-kW  total  power  at  10%,  20%,  40%,  60%,  and  100%  (of  total  slab  thickness)  liquid 
layer  retention  on  front  surface. 
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Appendix  A 

PROGRAM  DESCRIPTION 


The  program  is  set  up  so  that  with  an  initial  terminal  question-and-answer  period  the 
user  can  select  a material,  an  arbitrary  slab  geometry,  an  arbitrary  beam  profile,  and  an 
arbitrary  run  time.  An  internal  library  contains  all  the  material  parameters  and  allows  the 
user  simply  to  select  a materia]  without  the  aid  of  reference  tables.  The  target  thickness, 
length,, and  width  are  input  in  centimeters.  Beam  characteristics  are  input  as  total  power 
in  watts,  diameter  in  centimeters,  angle  of  incidence  in  degrees  from  normal,  and  selection 
of  Gaussian,  “top-hat”,  or  user-specified  experimental  profile.  The  ambient  or  starting 
temperature  is  input  in  degrees  centigrade,  and  boundary  conditions  of  heat  sinking  to 
ambient  or  insulating  among  edges  and  back  surface  are  selected.  Output  graphics  are 
set  up  by  inputting  total  run  time  and  printout  time  interval  for  slab  temperature  arrays 
and  plot  and  either  front  or  rear  surface  to  be  plotted. 

The  programing  is  done  in  Fortran  IV  for  improved  computation  speed  over  Basic. 
The  ratio  of  run  times  for  the  two  programing  languages  on  the  PDP-11  is  approximately 
40  for  a program  of  this  type.  The  entire  program  is  core  resident  on  a 28K-PDP  11/40, 
including  the  plotting  routine,  and  will  cycle  completely  through  one  time  step  (once 
through  program  loop)  in  1 s for  premelting  and  1.3  s for  postmelting. 
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PROGRAM  LISTING 


0001  DIMENSION  TO  ( 10 > 10.  1 1 » 2 ) » 1X0  (8.8  .9 ) .P<  10 » 10 ) . ,IH2(  10.  10) 

1 * I BUT"  < 1200)  .0(10) 

0002  10  FORMAT !'  1 ... ALUMINUM » 2 ...COPPER.  3 . . .MOLYBDENUM. ' > 

0003  20  FORMAT ( ' 4 ...IRON.  5 ...NICKEL.  6 ...’ITANIUM.  7 . . . AL < 2024  ) ' ) 

0004  30  FORMAT  < ' INPUT  * FOR  TARGET  MATERIAL') 

0005  40  FORMAT (14) 

0006  50  FORMAT ( ' MATERIAL  THICKNESS  (CM)') 

000/  60  FORMAT  (F'20.5) 

0008  70  FORMAT ( ' TARGET  LENGTH  (CM)') 

0009  80  FORMAT < ' TARGET  WIDTH  (CM)') 

0010  90  FORMAT ('  AMBIENT  TEMPERATURE  (O') 

0011  100  FORMAT < ' INCIDENT  POWER  (WATTS)') 

0012  110  FORMAT ( ' BEAM  DIAMETER  (CM)') 

0013  120  FORMAT ('  ANGLE  OF  INCIDENCE  (DEG)  FROM  NORMAL ' 1 

0014  125  FORMAT ( ' PREMELT  ABSORPTION  COEFFICIENT ' » 

0015  126  FORMAT ('  POST-MELTING  ABSORPTION  COEFFICIENT ' ) 

0016  130  FORMAT ('  TIME  OF'  RUN  (SEC)') 

0017  135  FORMAT ( ' PRINT  OUT  INTERVAL  (Si: O') 


0018  140  FORMAT < ' PRINT  X-AXIS  INTERCEPT  WITH  Y-Z  -LANE -1-9  > 

0019  145  F ORMAT  ( ' LIQUID  LAYER  RETAINED  (PERCENT  OF  ■ 4 1 CK NFS’S  .• 

0020  150  FORMAT!'  SELECT  MODE  #'} 

0021  151  FORMAT!'  0 ...FRONT  SURFACE  PLOT') 

0022  152  FORMAT ( ' 1 ...REAR  SURFACE  PLOT') 

0023  155  FORMAT!'  1 ...EDGES  INSULATED  .BACK  SURFACE  I.N  >ULA  i‘F.0  ' 

0024  160  FORMAT  ( ' 2 ...EDGES  HEAT  S INKED  1(3  AMECf.T-.rv  I » BACK  SURFACE 

1 INSULATED  ' ) 

0025  165  FORMA!  ( ' .5  ...EDGE.S  INSULA?  E.D  . BACK  SURFACE  He  fil  :>  i ME  " U ■ \ 

1 AMBIENT') 

0026  1.70  FORMAT  < ' 4 ...EDGES  AND  BACK  SURFACE  HEAT  51NKEB  Tf)  aHB-iEn 

0027  172  FORMAT ( ' INCIDENT  BEAM  PROFILE') 

0028  174  FORMAT ('  1 ...  GAUSS (AN ' > 

0029  176  FORMAT ( 2 ...TOP  HAT') 

0030  178  FORMAT! ' 3 .. .EXPERIMENTAL  PROFILE  ) 

0031  180  FORMAT < ' FRONT  SURFACE  ) 

0032  1 82  FORMAT ( 10 I 6 ) 

0033  184  FORMAT  ( ' EL  IF  MEN  ? ' . 1 4 . INTu  i WA  : t i)/(‘ri/CM ' 

0034  186  FORMAT < ' BURN  THROUGH  HAS  OCCURRED  A!  POSITION i 4. 14 > 

0035  IBS  FORMAT!'  AND  AT  TIME  ' . I-.I.0 . 5.  ' S..CONDS-  1 

0036  190  FORMAT! ' TIME  IN  SECONDS  " . F10.5) 

0037  191  FORhA  T ( ' FRONT  SURFACE') 

0038  1 92  FORMA T ( 1 0F8 . 2 ) 

0039  193  FORMAT!'  BACK  SURFACE  > 

0040  194  FORMAT!'  X-Z  SLICE  AT  Y-'.K!) 

0041  195  FORMAT!'  RUN  COMPLETED') 

0042  197  FORMA  I ( E 8 . 0 . F8 . 3 . FB , 5 . F8 .3) 

0043  TYPE  10 

0044  TYPE  20 

0045  TYPE  30 

0046  ACCEPT  40.111 

0047  GO  TO  (200. 210.220. 230.240.241. 242) .11 1 

0048  200  AO-660. 

0049  E 0= . 222 

0050  F2 - . 497 

0051  A 3=2 . 7 

0052  A4=94 . 

0053  E 1 = • 7 , 7 IE-5 

0054  F3--9 ,0E-4 
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ooss 

F 4=0  • 

0056 

F5=30 . 

0057 

F V = 2 * 4 3 

0058 

GO  10  300 

005V 

210 

00=1083. 

0060 

F0--.915 

0061 

F2= . 918 

0062 

03=8.89 

0063 

04-42. 

0064 

F 1=2 . 4E-5 

0065 

F3=-l .198-7 

0066 

F 4=0 . 

0067 

F3=18. 

0068 

F 9 = 8 . 2 1 7 

006V 

GO  TO  300 

0070 

220 

00-2610, 

0071 

F0= . 06162 

0071! 

F 2 = . 346 

00  73 

03=10.2 

00  74 

04-131 . 

0075 

Fl=2 , 2E-5 

0076 

F3-=-3 . 46E-‘ 

007  7 

F 4=0 . 

0078 

F5=10. 

007V 

F 9=8.1 

0080 

GO  ro  300 

0081 

230 

00=1535. 

0082 

F0= . 1060 

0083 

F2= . 108 

0084 

0 5=7 . 85 

0085 

04=65 , 

0086 

FI =9 . 6E-5 

0087 

F3=“l , 08E- 

0088 

F 4=0 , 

0089 

F5=l 8 . 

0090 

F 9=6 , 88 

009  L 

GO  ro  300 

0092 

240 

00= 1 453 ♦ 

0093 

F0=, 1095 

0094 

F2= . 1423 

0095 

03=8 . 75 

0096 

04=73 . 

0097 

FI =5 . 49E-5 

0098 

F3=-4 , 5 61:.— 

0099 

F4=0 . 

0100 

F5*0 . 

0101 

F 9 = 7 . 9 

0102 

GO  TO  300 

0103 

241 

00=1690. 

0104 

F0= . 1 39 

0105 

F2= . 0372 

0106 

03=4.54 

0107 

04-103.9 

0108 

F1=0 . 

0109 

F3—4.E-A 

0110 

F 4 = 0 , 

0111 

F5-21 . 

f 
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0112 

F9=4»09 

0113 

GO  TO  300 

0114 

242 

A0=630 . 

0115 

F0=  .215 

0116 

F2-.334 

0117 

03=2.77 

0118 

04=95.6 

0119 

Fl=7 , 71E-5 

0120 

F3=9.0E-4 

0121 

F4=0 . 

0122 

F5=30 , 

0123 

F9-2 . 43 

0124 

300 

TYPE  50 

0125 

ACCEPT  60  >06 

0126 

TYPE  70 

0127 

ACCEPT  60? 07 

0128 

TYPE  80 

0129 

ACCEPT  60? A8 

0130 

TYPE  90 

0131 

ACCEPT  60.A9 

0132 

TYPE  100 

0133 

ACCEPT  60? BO 

0134 

B0=B0/4 . 54 

0135 

TYPE  110 

0136 

ACCEPT  60  ? B 1 

0137 

TYPE  120 

0138 

ACCEPT  60? B? 

0139 

TYPE  125 

0140 

ACCEPT  60  ? A 5 

0141 

TYPE  26 

0142 

A C C E P T 6 0 9 A 

0143 

TYPE  130 

0144 

ACCEPT  60 ? B3 

0145 

TYPE  135 

0146 

ACCEPT  60  ? DO 

0147 

TYPE  140 

0148 

ACCEPT  40?. 1.0 2 

0149 

C9  = 0 . 

0150 

TYPE  145 

0151 

ACCEPT  40? I HI 

0152 

IH1= IH1/10 , + , 5 

0153 

TYPE  150 

0154 

TYPE  151 

0155 

TYPE  152 

0156 

ACCEPT  40  ? I PL 

0157 

IPL= 1+10+IPL 

0158 

TYPE  155 

0159 

TYPE  160 

0160 

TYPE  165 

016.1 

TYPE  170 

0162 

ACCEPT  40? IBS 

0163 

TYPE  172 

0.164 

TYPE  174 

0.165 

TYPE  176 

0166 

TYPE  178 

0167 

ACCEPT  40? IB6 

0168 

B7*AB/9 . 

8.T.  HANLEY 


0169 

89=81/87 

0190 

88=07/9. 

0171 

C0=Bl/B8 

0172 

GO  TO  <330 >320 >310) > 186 

0173 

310 

type:  180 

0174 

no  312  1=1 >10 

0175 

no  3.11  j=i.io 

0176 

IBUF < J)=J-1+10*(I-1 ) 

0177 

311 

CONTINUE 

0178 

TYPE  .1 82  > < 1 8 UP  < C. ) > L=  1 > 1 0 .) 

0179 

312 

CONTINUE 

0180 

80  314  1=1 >10 

0181 

DO  314  J=l>10 

0182 

R0= (1-5.5) **2/89/89+  < J-5 . 5 > **2/C0/C0 

0183 

IF  < RO  .GT.  ,25)  GO  TO  313 

0185 

K*  1'-- 1 + 1 0*  < J~  1 ) 

0186 

TYPE  184  > K 

0187 

ACCEPT  60>F6 

0188 

P<I> J)=F6/4.54 

0189 

60  TO  314 

0190 

313 

P(I>J)=0. 

0191 

314 

CONTINUE 

0192 

GO  TO  340 

0193 

320 

C4=4. *80/3. 14159/81/81 

0194 

DO  322  1=1 >10  1 

0195 

DO  322  J=l>10 

0196 

R0=< I -5. 5) **2/89/89+ ( J-5 ,5) **2/C0/C0 

0197 

IF  (RO  . GT , .25)  GO  TO  321 

0199 

P ( I > J ) =C4 

0200 

GO  TO  322 

0201 

321 

P(I>J)=0. 

0202 

322 

CONTINUE 

0203 

GO  TO  340 

0204 

330 

DO  332  1=1 >10 

0205 

DO  332  J=l>10 

0206 

R0=2 , * ( < < 1-5 ,5) *87 ) #*2+ ( ( J-5 . 5) *88 ) **2 ) ** . 5 

0207 

IF < RO  .GT.  81)  ‘GO  TO  331 

0209 

E0=-2.*R0/B1 

0210 

P ( I > J ) =4 , 287*80*EXP ( EO ) /B 1 /B1 

0211 

GO  TO  332 

0212 

331 

P(I>J)=0. 

0213 

332 

CONTINUE 

0214 

340 

m=no 

0215 

C1=0, 

0216 

DO  344  1=1 >10 

0217 

DO  344  J=l>10 

0218 

DO  344  K=1 > 1 1 

0219 

TO  < I > J>K> 1 )=A9 

0220 

T0< I > J >K>2 )=A9 

0221 

IF ( N > GT , 9 ) GO  TO  344 

0223 

I F ( I , GT . 8 . OR . J , GT  . 8 ) GO  TO  344 

0225 

IXO(I> J>K)=0 

0226 

344 

CONTINUE 

0227 

C2= > 2*  < FO+F 1 *A9 ) *A3/ < F2+F3* < A9-F5 ))/< 1 . /B8/B8+ 1 , /*7 787 
1 +100./A6/A6) 

0228 

86=02/87/87 

0229 

C0=C2*100 . /A6/A6 
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02  50 

C 3 C2/B8/P8 

0231 

C4--2,*C2*<  1 ./B8/B8+1 ./B7/B7+100./A6/ 

0232 

THET=3. 14159*82/180. 

0233 

63-COS  ( rut:  T ■ *A5*A6/5  . 

0234 

B4  -C2*10000,/A4/A3/B8/B« 

0235 

U C 2 * 1 0 0 0 0 r / 0 4 / A 3 / P 7 / B 7 

0236 

P a ■■■■  r 2 * 1 0 0 0 0 00  / 0 4 / 4 3 / a />  / a a 

0237 

BO -0  s 

0238 

DO  400  -6=1  . 10 

0239 

DO  400  1=1? 10 

0240 

IH2< I»  J)  = l 

0241 

DO-P  *:  I ? J ) # B7* B8*4  , 54+ BO 

0242 

400 

CONTINUE 

0243 

500 

ICt*8*Ci+2 

0244 

DO  505  K = 2 j .1 C 1 

0245 

DO  505  J-=2  ? 9 

0246 

DO  505  1=2 r 9 

0247 

IF<T0<I>J»K»1>.LE.A0)  GO  TO  505 

024V 

IFUXOa-l.J-l.K-l)  * NE  > 0)  GO  TO  50 

025.1 

TO  (I1J1 K • :!.  ) =A0 

0252 

1X0<  1-1 • J-l  ■ K - 1 > = i 

0253 

IF (Cl .GT. .5)60  TO  505 

0255 

Cl  = l . 

0256 

G 3 = G 3 # A M 5 / A 5 

0257 

A 5= AM 5 

0258 

505 

CONTINUE 

0259 

DO  510  J=.L  • 10 

0260 

DO  510  1 = 1 • 10 

0261 

I3=IH2  ( I y .J)  + :l 

0262 

10(1-  J • 1 3- :!.»!)  = < G3*P  ( I » J ) / < F2+F3*  ( TO 
1 4*  T O ( I • J ? 1 3 * 1 ) - TO  ( I j Jj  13+1 . 1 .1  >/3. 

0263 

510 

CONTINUE 

0264 

DO  525  K=2  > .10 

0265 

DO  525  J=2  !>  9 

0266 

DO  525  1=2 j 9 

0267 

IF  < TO  < I » J f K » 1 ) ♦ LE  » AO ) GO  TO  522 

0269 

F8=F9 

0270 

GO  TO  524 

0271 

522 

F8=A3 

0272 

524 

Al-  < F2+F3#  ( TO  CUJ.Kil)  -F5  > ) /FQ7  ( FO+F 

02/3 

to  -•  A1*D6*  ( TO ( I * J+l  • N r i ) +T0  < I !*  J-l « K • 1 

0274 

E.  I = A1#C0#  ( ! 0 < I ■ J »K+1  r 1 ) + TO  < 1 t J * K“  1 » 1 

0275 

E2=AL  *C3*  •:  TO  < I + i ? J « K ? 1 > + \ <>  • 1-1  ? J ? K ? 1 

0276 

T0(I» JfK» 2 ) = E0+E1+ E 2 + Cl . ~ A 1 * C 4 > * 1 0 < I 

0277 

525 

CONTINUE 

0278 

GO  TO  ( 530  ■ 535 , 540 , 550 ) j .1 B5 

0279 

530 

DO  532  J=1.10 

0280 

DO  532  1=1 ? 10 

0281 

l 0 • l j J • 1 1 • 2 > 1 0 ■ 1 • J ? lOr  2 > 

0282 

532 

CONI INUE 

02  £3 

540 

DO  3 K-2  ? 1 .1. 

0 2ci-» 

DO  533  J--1 . 10 

0285 

1 O'  1 » J ? K • J.  ' » I 0 < 2 • j . |\  , 2 :> 

0286 

i 0 1 0 j J j l\  .•  2 ' -•  f 0 < 9 j J * |\  » : : ) 

026/ 

533 

CON  T 1 N U C 

0288 

DO  534  l'  . . U 

026  v 

00  534  (=1 »10 

-F5) )+ 


< J » K » 1 ) ) 


...l  ? K * .t  i 
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i 0 ( I ? 1 ?K 
TO'.  I ? 10? 

i nJ  1 j.  HUt: 

GO  ro  55 


i .1  HUt: 

t r 1 C i : . ? t , 5 < GO  t ij  600 
.180  K=2?10 
HO  580  J - ■- 
DO  580  1=2? 9 
IF<IX0< 1-1 ? J-1?K~1 ) . L E , 
C7“F2+F 3* < TO < I ? J ? K ? 1 ) -F 
E0= TO  < I -1 ? J ? K ? 1 ) +T0 ( I + 1 
X XO  ( I -- 1 ? J- 1 ? K-l  > = I XO  < X - 


? J ? K ? :l  )-2*T0<I  ? JrK?  X ) 

1 < J-l ?K~1 M-C7*E0*B4+ . 5 


XI-  tlXO'I-l? J~1?K~ 
IX0<I~1?  J-l-K-l) 
XF  < ixo ( 1-1  ? J-l  ?!<:■- 
I. vox  X-:!.  ? J-l  ?K-1  )* 
XF ( K~IH2< I ? J ) . LT . 


X T , G T ? 0 ) GO  TC 

i 

1 ) iLE . 10000)  X 

-1 

IH1M30  TO  580 


X H2 ' X ? J ) = I H2 < X ? J > + X 
GO  TO  580 
T O ( I ? J ? K ? 2 ) =A0 

E0=1  0 ( I ? J-l  ? K ? 1 > +T0  ( I ? J+l  ? K ? 1 ) -2#T0 < I ? J ? K » 1 > 
IXOTI-l  ? J-l  ,K-1  )=I  XOU-l  , J-l  ,*K— 1 >+C7*E0*B5+.5 
IF ( X XO < I — 1 ? J— 1 » K— 1 ) * GT iO)  GO  TO  562 
IX0(I-1?J-1?K-1)=1 

IF  < IXO  < 1-1  ? J-l  ? K~1 ) ?l.E . 10000  ) GO  TO  564 

IXO( 1-1 ? J-l ?K-1 )=-l 

IF  < K-IH2  < I ? J ) . L.T  . IH1 ) GO  TO  580 

IH2<I?J)=IH2<I?J)+1 

GO  TO  580 

T0< I ? J?K?2)=A0 

EO=TO ( I ? J ? K-l ? 1 ) +T0  < I ? J * K+l ? 1 ) -2* T 0 ( I ? J ? K ? 1 > 
1X0<  X-l ? J-l ? K-l ) = IX0<I-1? J-l?K-l)+C7*E0*B6+.5 
IF(IX0<I-1?J-1?K-1) .GT.O)  GO  TO  572 
IX0(I-1? J-lfK-l)*l 

IF  < IXO  ( 1-1  ? J-l  ? K-l  ) . I..E  ? 10000 ) GO  TO  574 

IF <K.E«. 10)00  TO  573 

IX0(I-1 f J-l f K-l )*-l 

IF ( K-IH2  < I ? J ) . LT . IH1 ) GO  TO  580 

IH2(I? J)=IH2<I?J>+1 

GO  TO  580 

TYPE  186  ? I ? J 

TYPE  188?C9 

STOP 

T0<I?J?K?2) =A0 

CONTINUE 

C9=C9+C2 

IF  < C9 . L T . D1 ) GO  TO  700 

D1=D1+D0 

TYPE  190?C9-C2 

TYPE  191 

no  620  1=1? 10 

no  615  J=l?10 

IF(1.LT.IH2(I?J))G0  TO  612 
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0361 

01J)=T0<I.J.1.1> 

0362 

GO  TO  615 

0363 

612 

0< J)  = l .8+9 

0364 

615 

CONTINUE 

0365 

TYPE  1 92 » 1 0 1 L > . L. - 1*10) 

0366 

620 

CONT 1NUE 

0367 

TYPE  193 

036B 

80  630  1=1 , 10 

0369 

TYPE  192„<T0<X.L»11j1)»L=1.10) 

0370 

630 

CONT INUE 

0371 

TYPE  194.102 

0372 

DO  635  1=1 . 10 

0373 

80  631  K=2 . 1 1 

0374 

IF  < K . LT , I H2 1 1 » 182  > ) GO  TO  632 

0376 

0<K-1 >=T0< I . I82.K . 1 ) 

037  7 

GO  TO  631 

0378 

632 

0(K-1 )=1 ,E+9 

0379 

631 

CONTINUE 

0380 

TYPE  192. <0<L>?L= 1.10) 

0381 

635 

CONTINUE 

0382 

E0=-500 . 

0383 

E2=50000. 

0384 

80  640  J=l.l<) 

0385 

DO  640  1 = 1 . 10 

0386 

IF1T01I. J.IPL.l > ,LT.EO)GO  TO  636 

0388 

EO=TO( I » J» IPL  » 1 ) 

0389 

636 

IP(TO< I > J» IPL. 1 > . GT .E2)G0  TO  640 

0391 

E2=T0<  I » J»  IF'L.  1 > 

0392 

640 

CONTINUE 

0393 

IP<E0-E2.GT. .1)G0  TO  800 

0395 

El=,  1 

0396 

GO  TO  850 

0397 

800 

IF ( E0-E2  »GT . 1 . ) GO  TO  801 

0399 

El  = l. 

0400 

GO  TO  850 

0401 

801 

IF1E0-E2.GT ,5. )G0  TO  802 

0403 

El=5. 

0404 

GO  TO  850 

0405 

802 

IF1E0-E2  »GT , 10 , )G0  TO  803 

0407 

El =10 t 

0408 

GO  TO  850 

0409 

803 

IF (E0-E2  , G T . 50 , )G0  TO  804 

0411 

El =50 i 

0412 

GO  TO  850 

0413 

804 

I F < E0-E2 . GT , 1 00 , > GO  TO  805 

0415 

El =100 , 

0416 

GO  TO  850 

0417 

805 

IF ( E0-E2 , GT , 500 , ) GO  TO  806 

0419 

El =500, 

0420 

GO  TO  850 

0421 

806 

IF1E0-E2.GT, 1000, )G0  TO  807 

0423 

El =1000, 

0424 

GO  TO  850 

0425 

807 

El =3000, 

0426 

850 

SCA=125 , /El 

0427 

CALL  INITlIEiUF.il  20) 

0428 

CALL  SCAL10, ,0, .970, . 1000. ) 

21 


S.T.  HANLEY 


I 


0429  CALL  APNT  *.  ! 00  •.  - 1 00  . - i • - j -O  * 

04  <0  CALL  UEC  T < 00')  - 750 * * - i . a . - j w i . 

04  5 1 • ALL  APNI  < 1 00  • - 240  * ~ 3 - - ! - 0 » 

04  32  CAL  L VEC  I < 0 ; * - !.  4 0 * * - 1 * a • | « • v > 

043  3 C ALL  OLL  T < 500  • - Q * ? - 1 * a - - .1  - 1 » 

0434  DQ  64':.  I 1-5 

0435  CALL  APN  I < I * 1 00  » +50  : - 90  -a  • ••  • 

0436  CALL  UECT  < 0*10  - ~ 1 v 6 - - 1 - ' • 

0 4 3 7 6 4 5 C ONT IN U E 

0438  DO  a>47  .1.  • ! * 5 

0439  CALL  APNT <1*1 00 , +40 , - 1 * ! 50  +25 , - ~.L  * ~6  * 0 > 

0440  CALL  VECT/10.  *0  - -l*6*-lf  I.  > 

044.1  647  CONTINUE 

0442  DO  650  !~;l*6 

0443  CALL  APN  T < 90 » • 1 00  ♦ + » 1 - ! • * 25  * - - 1 * -“6  • 0 > 

0444  CALL  UEC  f < 1 0 « • 0 ♦ * - 3 <.6*1-  ! > 

0445  650  CONTINUE 

0446  DO  660  J • L * 10 

0447  C A L.  L APN  T ( 5 0 . + J *50#  - < 1 0 <!*.!•  I P L • 1 > - E 7 ' # 0 C A •*  ^ 5 t V J * '? 5 » - - .!.  ••  - ' • - - I ) 

0448  DO  460  Llf? 

0 449  CAL  L U E C 7 < 5 0 . • « T 0 < !..  + 1 - J • I p L - .1 .*  ••  T 0 ' L * J * I P > - .!  .*  > * S C A - - !••'••••  1 » '!  • 

0450  660  CONTINUE 

0451  DO  662  I-l-.IO 

0452  CALL  APNT  * 50  » + j *50  ••  < Vo  • I - .!.  - 1 pl  . 1 . -p2  *8C A+ 1. 00  • - - .»  ? - • - -•  .1  - j.  » 

045 3 DO  662  L ~ 1 * 9 

0 454  CALL  U E C T ( 5 0 * • <1  0 t * L + :l.  - 1 P L - 1 > - T 0 ( I -L-  j!  P L - ? • ■ * 8 C A + 7 5 5 - - 1 - i > 

0 4 5 5 662  C 0 N T I N U E 

0 4 5 6 C A L l A p NT*.  7 0 - 2 5 0 , * - 1 - - 6 • 0 > 

0457  CALL  TEXT < ' T+C) ' > 

0458  C A L L A p N T ( "3  0 0 • - 2 0 - - .!.  - - 6 * 0 > 

0459  1 1-  < I PL  * LT  *5 ) GO  TO  665 

*>461  C ALL  I E XT<  ' R E A R S U R cr  A C F C E N V j.  ft  hi' T E R S ' > 

0462  GO  TO  670 

0463  665  CALL  TEXT < ' FRONT  SURFACE  CENTIMETERS') 

0464  670  CALL  APNT <0 • *725  * *-l *-6 *0 > 

0465  CALL  TEX T < ' 3-D  TEMPERATURE  PROFILE  ' • 

0466  DO  675  I - J.  • 5 

0 4 6 7 CALL  APNT < I* 1 0 0 * - 1 5 * * 7 0 < - - 1 - - 6 - 0 » 

0468  CALL  NMBR  < I * 2*B8* < I - t 5 ) * " F7 < 2 ' > 

0469  675  CONTINUE 

0470  DO  680  1-2-5 

0471  CALL  APNT ( < I- , 5 > * 100 • • 1*150 * +15 < -~! *-6*0 > 

0 4 7 2 C A L L N M B R < I + 5 * 2 * B 7 * < I - ♦ 5 > - y F 7 1 2 ' > 

04  7 3 6 8 0 C 0 N T I N U E 

0474  DO  685  1-1*6 

0 4 7 5 C A L L A P N r < o * * 9 0 • + • I ~ 1 > * 2 5 < * - 1 - - 6 - 0 > 

04/6  J R < E 1 « G 7 » 5 0 , ) G 0 T n A P 4 

04 78  CALL  NMBR • I + I 1 - U - 1 > * 25  /SCA+E2 • ' F6  « 2 ' > 

0479  GO  TO  685 

0 480  684  C ALL  NMBR  *.  I + 11  * < I - 1 > * 2 5 * / S C A 4- E 2 - F6  * 0 ' « 

0481  685  CONTINUE 

0482  CALL.  APNT  < 0 . - 700  *-i*-6*0) 

0483  CALL  TEXT ' ' PREPARED  BY  DR  8 ♦ T ♦ HANLEY  N R « ") 

0484  CALL  APNT < 0 , • 675 ♦ * -1  * -6  * 0 ) 

0485  CALL  1 EXT (' MATERIAL ? ' > 

0486  CALL  APNT  < 1 4 1 » • 675  ; • - .!.  ? - 6 - 0 > 

0487  G 0 r 0 <690*691  -69 2 - 6 9 7 - 6 9 4 - 8 1 0 » 8 !. :!.  > * il  3.  .1 
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0488 

690 

CALL 

TEXT ( ' ALUMINUM'  > 

4 H?  ^ 

BO  T(] 

i 695 

04  VQ 

691 

CALL 

TEXT ( ' COPPER ' > 

0491 

GO  RJ 

i 695 

. j 4 V 2 

692 

CALL 

TEXT ( 7 MOL Y BO  ' > 

uJs>  A 

GO  TO 

i 695 

04  V 4 

693 

CALL 

TEXT ( 7 IRON ' > 

04  VS 

GO  TO 

i 695 

0496 

69  4 

CALL 

TEXT  < NICKEL'  > 

04  vy 

GO  TO 

| 695 

04^8 

810 

CALL 

TEXT  < 'TITANIUM' ) 

0499 

GO  TO 

i 695 

0500 

811 

CALL 

TEXT ( 7 AL  < 2024  > 7 ) 

0501 

695 

CAL  L 

AF’NT  ( 0 » f 650  • y-1 

0502 

CALL 

TEXT  < ' BEAM  POWER  5 ' ) 

0503 

CALL 

APNT  < 156 . r 650 , .-1*~6.0> 

0504 

CALL 

NMDR  < 21 y BO • ' F 8 . 0 ' ) 

0505 

CALL 

APNT ( 290  * 9 650  « • - 1 • -6 • 0 ) 

05 06 

CALL 

TEXT ( 'WATTS' ) 

0507 

CALL 

APN ! (0*?  625 » f - 1 • -6  r 0 > 

0508 

CALL 

TEXT < ' IRRADIATION  TIME  t ' > 

0509 

CALL 

APN  T < 223  * * 625  « • -1 9 -6 • 0 ) 

0510 

CALL 

NMBR ( 22 9 C9-C2 • 7 F8 • 3 ' > 

0511 

CALL 

APNT  (350,  *625.  *-l*-6*0> 

0512 

CALL 

TEX  T < ' S ' ) 

0513 

CALL 

APNT  < 0 . .600, .-I .-6.0) 

0514 

CALL 

TEXT  < ' ABSORPTION  COEF  , ' ) 

0515 

CALL 

APNT ( 238 , i 600 , .-I .-6.0) 

0516 

CALL 

NMBR ( 23 * A5 * 'F8.5' > 

0517 

CALL 

APNT  (0,  ..575,  *~1  .-6.0) 

0518 

CALL 

TEXT ( ' BEAM  DIAMETER? ' ) 

0519 

CALL 

APNT (193. r575, r-1 .-6.0) 

0520 

CALL 

NMBR(24.B1 * 'F8,3' ) 

0521 

CALL 

APNT ( 335 . ? 575 , * - 1 * -6 * 0 ) 

0522 

CALL 

TEXT ( 'CM' ) 

0523 

700 

IF  ( C9 

’ . GT , B3 , OR  > C9 , EQ , B3  > GO  TO 

0525 

BO  7.10  K-l  r 1 1 

0526 

DO  710  J—l f 10 

0527 

DO  710  1=1 » 10 

0528 

TO  ( [ * 

J . K f 1 ) -TO  ( I * J t K * 2 > 

0529 

710 

CONTINUE 

0530 

GO  TO 

I 500 

0531 

750 

TYPE 

195 

0532 

END 

